function [I,J,K] = indices(d)
%        [I,J,K] = indices(d)
m = (d+1)*(d+2)/2;
I = zeros(m,1);
J = I;
K = I;
idx = 1;
for i = d:(-1):0
    for j = (d-i):(-1):0
        I(idx) = i;
        J(idx) = j;
        K(idx) = d - i - j;
        idx = idx + 1;
    end
end